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Abstract 

The joint management of heat and power systems is believed to be key to 
the integration of renewables into energy systems with a large penetration 
of district heating. Determining the day-ahead unit commitment and pro¬ 
duction schedules for these systems is an optimization problem subject to 
uncertainty stemming from the unpredictability of demand and prices for 
heat and electricity. Furthermore, owing to the dynamic features of pro¬ 
duction and heat storage units as well as to the length and granularity of 
the optimization horizon (e.g., one whole day with hourly resolution), this 
problem is in essence a multi-stage one. We propose a formulation based on 
robust optimization where recourse decisions are approximated as linear or 
piecewise-linear functions of the uncertain parameters. This approach allows 
for a rigorous modeling of the uncertainty in multi-stage decision-making 
without compromising computational tractability. We perform an extensive 
numerical study based on data from the Copenhagen area in Denmark, which 
highlights important features of the proposed model. Firstly, we illustrate 
commitment and dispatch choices that increase conservativeness in the robust 
optimization approach. Secondly, we appraise the gain obtained by switch¬ 
ing from linear to piecewise-linear decision rules within robust optimization. 
Furthermore, we give directions for selecting the parameters dehning the un¬ 
certainty set (size, budget) and assess the resulting trade-off between average 
proht and conservativeness of the solution. Finally, we perform a thorough 
comparison with competing models based on deterministic optimization and 
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stochastic programming. 
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Nomenclature 

Sets 

t time index; 
k unit index; 
i storage index; 

Parameters 

minimum/maximum heat output of unit k] 
vJVk minimum/maximum power output of unit k\ 
rj^jrk downward/upward ramping limit for fuel input of unit k\ 

TkK minimum up-/down-time for unit fc; 

Ui/ui minimum extraction/maximum injection of heat from/to storage i; 
Sj/sj minimum/maximum heat content of storage d 
hxed cost for unit k when on; 

Cfc per unit fuel cost of unit k; 

start-up/shut-down cost for unit k; 
dt heat demand at time t. 

Decision variables 

Vkt on/off status of unit k at time t (binary); 
ykt start-up status of unit k at time t (binary); 

Zkt shut-down status of unit k at time t (binary); 

Qkt heat production from unit k at time t; 

Pkt power production from unit k at time f; 
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Uit heat injection into storage i at time 
Sit heat content in storage i at time t; 


1. Introduction 


Combined Heat-and-Power (CHP) plants are highly efficient units where 
exhaust gases from the electricity production process are used to produce 
heat (co-generation). That way, energy that would otherwise be spilled in the 
atmosphere is further exploited for residential or industrial heating purposes. 
It is estimated that CHP units based on the Combined Cycle or Open Cycle 
Gas Turbine (CCGT, OCGT) technology can reach an efficiency of 90%, 
while typical val ues are 60% for standard powe r-only plants based on the 
same technology ( Veerapen and Beerepootl. 201 ih . 

Another important feature of CHP plants is their flexibility. Gas-hred 
GHP units, just like their power-only counterparts, are able to change the 
level of their output (ramp-up/-down) relatively quickly. Furthermore, they 
enjoy relatively short start-up and shut-down time. Moreover, GHP units 
have two additional levels of flexibility compared to power plants. Firstly, 
they are typically coupled with heat storage facilities, i.e., insulated tanks 
where hot water can be stored with negligible losses. This allows for a shift 
in time between production and consumption of heat. Secondly, the ratio 
between the heat and power output is variable for the so-called extraction 
GHP plants. This gives them extra flexibility, in that the heat and power 
outputs can adapt to the demand and price of these commodities. 

As a result of their attractive properties in terms of efficiency and flex¬ 
ibility, co-generation is now experiencing a revival in view of the increasing 
commitment of industrialized countries to curb GO 2 emissions. Indeed, GHP 
units provide an efficient way of producing heat that c an be made essentially 


carbo n neutral by switching from fossil fuel to biomass fiFernandez Pales et al. 
2 OI 4 J) . Furthermore, owing to their flexibility, GHP units can support the 


integration of large-scale electricity production from variable and partly- 
predictable renewable sources, e.g., wind and sol ar, by providing back-up 


powe r when the output of the latter is insufficient fIVeerapen and Beerepoot 


20111 ). To unlock this potential, it is critical that the heat infrastructure and 
production assets are thought, planned and managed effectively and jointly 
with the other energy carriers (electricity, gas), the transport sy stem and the 
industrial assets in the same urban area ( Meibom et ahl . 2013 1. This work 
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addresses a fundamental part of this complex problem: the joint optimization 
of short-term operation of heat and power generation units. 

As far as electricity is concerned, our assumption is that the owner of 
the production assets elects the day-ahead market as the preferred floor for 
trading their electricity output. This is consistent with the current situation 
in many regional power markets, where t he largest pa. rt of the exchange of 
electricity takes place one day in advance fjWeberl. I 2 OIOI) . As the gate closure 
for day-ahead markets is typically set around noon on the day before delivery 
of electricity, participants must schedule electricity production, including the 
on/off status of the production plants (unit commitment). Similar arrange¬ 
ments requiring the day-ahead scheduling of heat production may also be in 
place. Indeed, the compatibility of the heat production plan with respect to 
the transmission and distribution constraints must be verihed in advance by 
the district heating network operator. 

As a result of the arrangements described above, the owners of CHP plants 
must determine the commitment of their units and the production schedule 
for heat and power with a signihcant advance in time. Many important 
parameters of the problem, however, are not known at the day-ahead stage. 
For example, this is the case of heat demand and electricity prices, both of 
which are only partly predictable. Another feature of the problem at hand is 
its multi-stage nature, which results from the fact that the unit commitment 
and production schedules to be determined span a whole day (typically with a 
granularity of one hour). This is particularly important as the studied system 
is inherently dynamic, owing to the presence of intertemporal constraints such 
as ramping limits and to the state-space representation of storages. 

In view of the stochastic and dynamic nature of the problem, the frame- 
wo rk of Affinely Adjust able Robust Optimization (AARO) is appropriate. 


see iBen-Tal et ahl (120091 ). In our view, this is one of the few mathematical 
programming techniques that reconcile an appropriate representation of the 
uncertainty with a rigourous modeling of the multi-stage nature of the system 
that does not violate the non-anticipativity of the sequence of decisions. 

The existing literature on the short-term operation of heat and power 
systems focused on the use of different optimization techniques. The prob¬ 
lem of optimizing the opera tion of a storage-CHP un it system is tackled via 
stochastic programmin g bv iPalsson and Raynl fll994l j and via stochastic dy¬ 
namic programming by iRavn and RvgaardI fjl994l ) . A deterministic approach 
for sch eduling a district heating system is studied bv lAringhieri and Malucelli 
(120031) . resulting in a Mixed Integer Linear Programming formulation of 
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the problem. O nly more recently did electricity markets enter the picture. 
Rolfsman f 2004f) used a stochastic programming approach to the problem 


including a model of the uncertainty in heat demand and day-ahead power 
prices. The stochastic programming a pproach has been enriched by fur- 
ther considering the intra-day stage by iDe Ridder and ClaessensI ( 20141) or 
by allowing for the subrnission of offering curves in the day-ahead market by 
Dimonlkas and Amelin ( 20141) . 


Recently, a number of papers studied various applications of robust opti¬ 
mization to power systems. The unit commitment problem has been studied 
via (two-stage) Adaptive Robust Op timization (ARO) consid ering uncer 


tainty in load (IBertsimas et al.L l2012h. win d power production Ijiang et ah 
(120121) or contingencies IStreet et ahl (1201 if) . The dispatch problem, where 
unit commitment variables are not consider ed, for both energy and r eserve 
has been studied with a similar approach by Zugno and Coneia ( 2013 ) with 
a vie w to immunizing the sol ution to uncertainty in wind power produc¬ 
tion. Warrington et ah ( 2013[) incorporated linear decision rules within a 
stoch astic optimization approach to a similar dispatch problem. Further¬ 


more, iLima et al.l (120151) proposed an ARO model for the short-term opti¬ 


mization of a virtual power plant under uncertainty in prices and wind power 
production. 

In view of the state-of-the-art described above, the contributions of this 
paper are the following. Firstly, we propose an Affinely Adjustable Robust 
Optimization (AARO) model for the unit commitment and dispatch problem 
for units participating in both the electricity and heat markets. Secondly, we 
set up a case-study using realistic data for unit parameters, heat consumption 
and prices obtained for the Copenhagen area, which accounts for ro ughly 20% 
of the total consumption of heat in Denmark ( Varmelast.dk . 20141) . Thirdly, 
we conduct an extensive numerical study where we assess the use of AARO 
for the problem at hand from the following perspectives: performance of 
purely linear decision rules against piecewise linear decision rules, sensitivity 
of results to changes in the parameters dehning the uncertainty set, and 
comparison with the competing approaches of deterministic optimization and 
stochastic programming. 

The paper is structured as follows. Section |2] introduces the deterministic 
version of the optimization model for unit commitment and dispatch of the 
heat and power system. Uncertainty in electricity prices and heat demand 
is then incorporated by reformulating the model as an AARO problem in 
Section [31 The numerical study based on data from the Copenhagen area is 
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presented and its results are discussed in Section |U Finally, conclusions are 
drawn in Section 13 

2. Deterministic modeling of heat and power system 

This section introduces the deterministic optimization model for the unit 
commitment and dispatch of a heat and power system. The first part of 
the section is dedicated to modeling of production units: back-pressure CHP 
units are dealt with in Section 12.11 extraction CHP plants in Section 12.21 
and heat-only plants in Section 12.31 Then, heat storages are described in 
Section ITTI Finally, Section 1231 loins the unit models in a common market 
framework. 

2.1. Back-pressure units 

Back-pressure units are a particular type of CHP units that can only 
produce heat and power at a fixed ratio. Let /C® be the set of units k of the 
back-pressure type, and t a time period within the optimization horizon. A 
back-pressure unit is characterized by a heat-to-power ratio r\ such that: 

Pkt = rl-qkt, yke}C^,t. ( 1 ) 

The heat output for the unit is limited both from below and from above 
by and respectively, via the following constraints: 

q^Vkt < qkt < qkVkt, yk eJC^,t. ( 2 ) 

The binary variable Vkt G {0,1} denotes the on/off status of unit k, forcing 
a null heat output when the unit is off. 

Given that the power and heat outputs are linearly dependent, there is no 
need to enforce upper and lower bounds for the power output, which follow 
trivially from ([T]): 

Pk = ^k-Qk' Vfc e /C®, (3) 

Pk = rl-qk, ykeJC^ . (4) 

We approximate fuel consumption, fkt, as a linear function of the heat 
and power output: 

fkt = plPkt + Plqkt yk , (5) 
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where and represent the fuel consumption per unit of power and heat 
output, respectively. 

Technical limitations restrain the change in operating regime (ramping) 
of a unit between consecutive time periods. The following constraints impose 
lower and upper ramping limits, and r^, on fuel consumption of unit h. 

-Lk < fkt - <fk Wk . (6) 

As an alternative, similar constraints could be imposed on the individual 
heat and power outputs rather than on the fuel consumption. 

The additional binary variables ykt,Zkt G {0,1} indicate, if equal to 1, 
that unit k is being started up or shut down at time t, respectively. They are 
related mutually and to the unit status, Vkt, via the following constraints: 

Vkt - Vk{t-i) - Vkt + Zkt = 0 , Wk , (7a) 

Vkt + Zkt<l , \/keK?,t . (7b) 

Constraint (l7b|) imposes that a unit cannot be started up or shut down at 
the same time. Note that such a constraint is not needed if the start-up and 
the shut-down operations have a strictly positive cost. 

Minimum up-time is imposed via the following: 




nUO 

'-k 


Vkt = Tl 

t=l 

TUt(t) 

^ ^ {^kr l/fct) ^ 0 


\^k e /C^ 


VA: e 


(8a) 


(8b) 


T = t 


The initial up-time is enforced by flSal) . Similarly, fl8bl) enforces min¬ 
imum up-time during the other time periods. Note that in the latter, we 
define 


rUf(t) =min{f + TU-l,r} , WkE)C^,t 


(9) 


where T indicates the final time period in the horizon and the minimum 
up-time for unit k. The constraints enforcing minimum down-time are spec¬ 
ular to (IHll and skipped here for the sake of conciseness. We refer the reader 


to 


Morales et al.l (120141) for further details. 


Finally, we define the cost of operating the back-pressure unit k at time 
t as the sum of fuel, no-load, start-up and shut-down costs: 

costn = Ckfkt + clvkt + cf^Vkt + cf^Zkt , VA: e /C®, t . (10) 
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2.2. Extraction units 

Extraction units allow for more flexibility in the ratio between heat and 
power outputs than back-pressure units. Let us consider a unit k of the set 
k G /C^ of extraction units in the system. 

The dehnition of the feasible region in the heat-power space is based on 
the following relaxation of the back-pressure constraint ([T]): 

Pkt>rl-qkt V/ce/C^,f. (11) 


Like in the case of a back-pressure unit, lower and upper bounds for heat 
production are dehned by ([2]). 

Fuel consumption is given by ([5]). For an extraction unit, we dehne upper 
and lower bounds on fuel consumption; 


l^Vkt < fkt < fk'^kt , Wk eK:^,t . (12) 

Note that, as a result of dehnition ([5]), flT^ implies that the feasible set in the 
heat-power space for each unit k G /C^ and at each time period t is included 
in the region between two parallel lines with the following slope: 

< = -4 , (13) 


When fuel consumption is at the upper bound, this slope represents the 
decrease in electricity output needed to produce an extra unit of heat without 
exceeding the upper bound in flT^ . H ence it i s often referred to as marginal 
electricity loss for heat production, see Weberl ( 20051 ). 

Constraints (I6])-(IH]), i.e., ramping limits, binary variable constraints, min¬ 
imum up- and down-time as dehned for the back-pressure unit in Section 12.11 
hold for an extraction unit too. Similarly, dehnition flT0|) of operating cost 
also holds. 


2.3. Heat-only production units 

Let us dehne unit k of the set as a unit only capable of producing 
heat. We impose lower and upper bounds on heat production ([2]), relations 
between binary variables (|7]), minimum up-time (IH])-([9]) as well as minimum 
down-time. 

Fuel consumption for these units is solely determined by heat production: 

fkt = pIqu , Vfc G /C“, t , (14) 
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Fuel-to-heat coefficient can in this case be interpreted as the inverse of 
the unit efficiency. 

Ramping limits on fuel consumption ([6]) still hold. In this case, they can 
be directly translated to ramping limits on heat production. Furthermore, 
the definition of the operating cost flTOj) is also valid for this type of unit. 

2.4- Heat storage units 

Let us consider a heat storage i from the set X of accumulator tanks 
installed in the system. 

The level of heat stored in the tank, su, is governed by the following 
state-update equation: 


Sit — Sip-1) + Uit , Vi E X,t , (15) 

where Uu represents the flow of heat from the production units into the tank 
at time t. 

We impose the following lower and upper bounds on the level of heat in 
the storage: 


Si < Sit < Si , Vi eX,t . (16) 

Typically, the lower bound is set to 0, while the upper bound is the thermal 
capacity of the storage. 

Moreover, we assume that the flow of heat to/from the storage tank is 
bounded from below and above: 

Ui <Uit<Ui, Vi El,t . (17) 

Finally, we enforce the condition that the final heat content in the storage 
be equal to the initial one through: 

SiT = Sio , Vi . (18) 

Note that this condition is needed to avoid that the storage be completely 
emptied throughout the optimization horizon. An alternative solution would 
consist in assigning an economic value to the heat stored in the tank at the 
last time period. However, we prefer the solution with terminal condition 
ffTSj) for two reasons. Firstly, the value of stored heat is not easily quantifi¬ 
able. Secondly, ffTSj) makes the comparison between competing optimization 
methods more transparent, as it ensures that the same amount of heat is 
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stored at the end of the horizon regardless of the modeling choices. In the 
alternative solution where a value is assigned to the heat stored at the end 
of the horizon, that comparison would be dependent on the arbitrary choice 
of the economic value for the terminal heat level. 


2.5. Heat and power market framework 

We consider an electricity pool with tii ne-varying prices set by the in¬ 
terception between supply and demand, see Morales et ah f 2014|) . As far as 
heat is concerned, we consider a market with a hxed price set by a regulator. 

We make the assumption that the owner of the CHP system is responsible 
for satisfying the demand of heat in a given urban area. This is enforced 
mathematically via the heat balance equation: 


Uit = dt , Vf , (19) 

keK iei 


where dt is the heat load at time t and K. = /C® U U /C^. 

The objective function is given by the market revenues from electricity 
sales minus the total costs for operating the units: 


P = {\Pkt - Cfc [pIpm + pWi] - clvkt - cf^Vkt - cf^Zkt] , (20) 

t,k^K 


where Xt is the electricity price at time t. Note that sales of heat are not 
included in the objective function, because the price for heat is hxed and so is 
the production over the optimization horizon, as a result of the combination 
of (1T8|) and flT^ . Consequently, revenues from heat sales are a constant, and 
hence they can be omitted from objective function (|20|) . 


3. Decision-making via Affinely Adjustable Robust Optimization 

The optimization model presented in the previous section is determinis¬ 
tic, in that it relies on the assumption that all parameters be known with 
certainty at the time of decision-making. However, this assumption is unre¬ 
alistic for the problem at hand, as unit-commitment and dispatch decisions 
must be made before the gate closure of the day-ahead markets described in 
Section [H i.e., when heat demand and power prices are unknown. In this 
section, we cast the problem as a two-stage optimization problem within the 
framework of Affinely Adjustable Robust Optimization (AARO). 
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3.1. Formulation as an Affinely Adjustable Robust Optimization problem 
In order to formulate the problem within the AARO framework, we shall 
hrst endow the decision-maker with the capability to perform recourse deci¬ 
sions, i.e., adjustable to the realization of the uncertain parameters. Let us 
collect all the hrst-stage variables in a vector x, and dehne the vector of 
recourse decisions as follows: 


V 

y 

z 


X = 


p , 

q 


s 

u 


ya 


Ap5 

Aq5 

Ass 

Au5 


( 21 ) 


Unit commitment variables belong to the set of hrst-stage decisions, as unit 
status cannot be changed on a real-time basis. The vector x further includes 
the day-ahead dispatch vectors p, q, s, u. Real-time adjustments with respect 
to the dispatch are included in the vector of recourse variables y^. 

With the dehnitions above, we can formulate a stochastic version of the 
optimization problem described in Section [2] as follows: 


Min. Ea {cjx-h gjya} 


x,y5 

s.t. Ax > b 


T^ + W^ys>K- 


(22a) 

(22b) 

(22c) 

(22d) 


We grouped all the constraints involving only hrst-stage variables into fl22bj) . 
while the ones including recourse variables are divided into equalities fl22cl) 
and inequalities fl22dD . These last two types of constraints are subject to 
right-hand side uncertainty stemming from the stochastic nature of the heat 
demand in (flU]) . The cost coefficients in the objective function in fl22al) are 
also uncertain as they include electricity market prices, see (I2U]) . We assume 
the same marginal cost for dispatch x and adjustments y^. However, the 
power dispatch is sold at the day-ahead market price, while the real-time 
adjustment is sold or bought at the balancing market price. We assume 
price consistency, i.e., that the e xpectation of t he ba lancing market price be 
equal to the day-ahead price, see Zavala et ahl ( 2014 ). 
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In order to put us in the framework of AARO, we make the following 
assumptions: 


Al. We require fl22c 


flBertsimas et ah 


l -fl22dl) be valid V(5 G W, where the uncertainty set 


20111 ) W is a bounded polyhedron described by the 


A2 


set of I linear inequalities L5 > 1. 

The uncertain parameters g^, h^, depend linearly on a random vec¬ 
tor (5, i.e., gs = g-|-G(5, Note that there 

is no loss of generality here, as we could simply redefine the uncertainty 
vector (5 as the concatenation of all the stochastic parameters in the 
model, and also the associated probability space accordingly. 

A3. The recourse decision is restricted to be an affine function of the 
uncertainty, i.e., = Y(5. 


After replacing the affine dependencies to the uncertainty in 
obtain the following: 

Min. E 5 {cjx + g^Y5 + (5^G^Y(5} 


s.t. Ax > b , 

T"x + W^^YS = h" + , 

mjn { (W'Y - H‘) > R - T‘x , 


yseu, 


we 

(23a) 

(23b) 

(23c) 

(23d) 


where fl23d|) is equivalent to requiring that (I22d|) be valid \/d eU (note that 
the min operator acts row-wise). With some further reformulations we can 
cast the problem as the following Mixed-Integer Linear Problem (MILP): 


Mrn^E^ {cj} X + g^YE^ {6} + tr {G’^Y (S 5 + E{(5}E{5}^) } (24a) 


s.t. Ax > b , (24b) 

T^x = R , (24c) 

W®Y = , (24d) 

A^l > R - T'x , (24e) 

L^A = (W'Y-H')^ , (24f) 

A G Mir > (24g) 

Y G M”^", r = 0 , V(z,i) G Z . (24h) 
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In order to obtain fl24ap . we exploited the invariance of the trace operator to 
cyclic per mntations of the operand s. The pairs fl24cl) - (l24dD are equivalent to 


see 


Warrington et al.l 


peranas 

I 


Note that fl24cp dehnes the day-ahead 
dispatch as the response to the nominal value of the right-hand side h®. Con¬ 
straints fl24el) - p^ are equivalent to fl23dP . They are obtained by applying 
strong linear duality row-wise to each of the constraints in fl23dp . as com¬ 
monly done in robust optimizat ion to determine the ro bust counterpart of a 
single inequality constraint, see iBertsimas et al.l (120111) . The (non-negative) 
dual variables dehned in each of these operations are collected into the ixm 
matrix A, where i is the number of constraints dehning the uncertainty set 
U and m the number of inequality constraints in fl23dp . Constraints fl24hp 
dehne some elements of the matrix Y be equal to zero. As explained in Sec¬ 
tion [331 this is necessary to guarantee the non-anticipativity of the solution, 
as well as to implement either purely linear or piecewise-linear decision rules. 

Note that the objective function in fl24al) includes the conditional mean 
of the day-ahead power price (included in c^), as well as the mean, ]E{<5}, 
and the variance-covariance matrix, S 5 , of the uncertainty 5. 

3.2. Definition of budget uncertainty set 

In this section, we briefly introduce the so-called budget uncertainty set, 
which is used throughout this paper to model lA in the AAR O model of the 
previous section. We refer the reader to iBen-Tal et al.l (120091) for a thorough 
introduction to the concept of uncertainty set. 

In this work, we consider budget uncertainty sets dehned by the following 
linear constraints: 


-1<S <1 , (25a) 

ll^lli = 1^1 < r . (25b) 

Inequalities fl25al) limit the uncertainty in a symmetric interval centered 
about zero and of radius equal to one. Note that the affine dependence be¬ 
tween the uncertain parameters in optimization model fl22D and the variables 
<5 allows us to consider any nominal values (g, h®, R) and interval sizes 
(which are ultimately set by the matrices G, H®, H^) for electricity prices 
and heat demand. Constraint fl25bD dehnes the budget of uncertainty for d, 
enforcing that its total absolute deviation from zero (i.e., summed across all 
the elements of the vector) is no larger than T. The latter is a parameter 
that can be tuned by the decision-maker. Larger values of T imply a more 
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robust solution, as it must be able to cope with uncertainty taking values 
within a larger set. Note that through the affine dependence of the parame¬ 
ters from d, fl25bD allows us to enforce rules such as: “during the considered 
optimization horizon, an uncertain parameter can have a deviation from the 
nominal value equal to the interval size during at most F hourly periods”. 

Because of the presence of the absolute value operator, constraint fl25bj) 
is nonlinear. However, it can be represented linearly by extending the uncer¬ 
tainty space with the inclusion of positive and negative parts of d: 


6' = 

'6' 

<5+ 

_ 

6 

max(<5, 0) 


d\ 


— min(5, 0) 


where max and min operate row-wise. A linear representation of the uncer¬ 
tainty set can be then obtained by replacing fl2^ with the following linear 
constraints: 


s = d+ -6- , (27a) 

0<(5+,5~<l, (27b) 

1^ (<5+ + (5”) < r . (27c) 


3.3. Definition of decision rules 

Decision rules are dehned by the matrix Y, which is a decision variable in 
optimization model fIM)) . Matrix Y dehnes the real-time adjustment of units 
and storage operation as a function of the realization of the uncertainty 5. In 
principle, the dehnition of Y in fIM)) may establish links between adjustments 
of operation variables and uncertainties whose value has not been revealed 
yet, hence destroying the non-anticipativity of the decision structure. Thus, 
additional constraints must be imposed to enforce that the elements of Y 
linking the readjustment to uncertainty unfolding at future time periods be 
zero. If there is one element of 6 per time period and the vector is ordered 
chronologically, this corresponds to requiring that Y is the concatenation of 
lower triangular matrices (one matrix per type of readjustment variable). 

Moreover, let us consider the extended dehnition 0261) of the uncertainty 
space including positive and negative variables. Similarly, we can dehne 
extended decision rules Y' including, besides the Y coefficients for 6, also 
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Y+ and Y” for 5’’" and S~, respectively. We can define purely linear decision 
rules by enforcing Y”*", Y“ = 0. Similarly, piecewise-linear decision rules are 
defined by setting Y = 0. 

Constraint fl24hD enforces that some elements of the Y matrix be zero 
under a choice of the index set 2 that is consistent with the two observa¬ 
tions above, i.e., it complies with the non-anticipativity of the solution and 
properly implements either linear or piecewise-linear decision rules. 

We conclude this section with an important remark. When using piece- 
wise linear decision rules, after replacing Y with Y' and <5 with S' in the 
objective function in ([24a]), we obtain (omitting zero terms): 

E, {cj} X + g^Y+E, [6+] + g^Y^E^ {6-} 

+ tr{G’^Y+E4(5+5^} + G^Y"E45^(5^}} . (28) 

Hence, the use of piecewise linear decision rules calls for a richer model 
of the uncertainty than with purely linear decision rules, as it requires the 
estimation of means and expectations of products between <5, <5'*' and <5“. 


4. Case study 


In this section, we present and discuss results obtained from an extensive 
numerical study assessing various features of the proposed model. 

We consider a system comprising an extraction and a back-pressure CHP 
unit, a peaker producing only heat and a heat storage. The technical pa¬ 
rameters for the production units are listed in Table [U Parameters for the 
CHP units are based on values for exi s ting u nits of the Copenhagen area 
as reported in IVarmeplan HovedstadenI (120141) . Notably, the back-pressure 
unit is the cheapest. However, it is less flexible than the extraction unit, 
which, besides having larger ramping limits and shorter minimum up- and 
down-time, can produce heat and power at different ratios. The heat-only 
unit is meant to be employed only as a back-up, i.e., during periods of peak 
demand. Hence, we assumed larger costs for operating this unit. This unit is 
also very flexible in terms of minimum and maximum heat output, ramping 
and minimum up- and down-times. We make the assumption that the back¬ 
pressure unit cannot modify the production schedule at the real-time stage, 
hence it is marked as “non-flexible” in Table [U This models a unit that for 
technical reasons (e.g., the type of fuel used) is too slow to perform changes 
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Table 1: Parameters for the units in the illustrative example 


Description and symbol 


Unit 

CHP 

extr. back-pr. 

Heat-only 

Power/heat ratio 

n 

- 

0.64 

0.28 

0 

Power loss for heat prod. 

Tv 

~ 

-0.13 

- 

- 

Fuel per electricity unit 

(pp 

- 

2.40 

2.40 

- 

Fuel per heat unit 


- 

0.31 

0.36 

1.09 

Min. fuel input 

f 

MWh 

120 

72.24 

0 

Max. fuel input 

f 

MWh 

631.20 

516 

1086.96 

Ramp-up limit 

r 

MWh/h 

150 

50 

1086.96 

Ramp-down limit 

L 

MWh/h 

150 

50 

1086.96 

Min. heat output 

Q 

MWh 

0 

70 

0 

Max. heat output 

q 

MWh 

331 

500 

1000 

Min. power output 

P 

MWh 

41.56 

19.60 

0 

Max. power output 

P 

MWh 

263 

140 

0 

Fuel cost 

c 

€/MWh 

24.16 

12.75 

93.96 

No-load cost 

cO 

€ 

0 

0 

2684.56 

Start-up cost 

cSU 

€ 

7382.55 

6 040.27 

0 

Shut-down cost 

^SD 

€ 

7382.55 

6 040.27 

0 

Min. up-time 


h 

2 

5 

0 

Min. down-time 

rj-^D 

h 

2 

5 

0 

Flexible 

- 

- 

yes 

no 

yes 


in its output level with short notice. We consider a storage tank with a heat 
capacity s = 2000 MWh and maximum in-/out-flow u = —u = 300 MW. 

We employ real-world data for electricity prices and heat demand in the 
Copenhagen region in Denmark. Power prices for Eas tern Denmark (DK 2 
area of Nord Pool) are publicly available for download at lEnerginet.dkl (120151) . 
We use heat consumption d ata from the western Copenhagen area, VEKS, 
available at Madsenl ( 2015 1. Note that the latter data refer to the time 
period between July 1995 and June 1996. Since the power market was not 
liberalized at that point in time, we paired heat load data with power price 
data from the corresponding week in the period 2013-2014. We rescaled the 
heat consumption data linearly so that the heat-only plant covers roughly 
0.5% of the total load over the year. We consider four weeks representative 
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of summer (the 4th week after the start of the dataset), autumn (17th), 
winter (30th) and spring (43rd). The winter week includes the largest hourly 
realization of heat load over the entire year. 

In the remainder of this section, we consider the following methods: 

1. Deterministic Optimization (DET): corresponding to the model de¬ 
scribed in Section [2] with expected values in place of Xt and dt] 

2. Robust Optimization with Linear Decision Rules (RO-LDR): imple¬ 
menting (124)) using purely linear decision rules; 

3. Robust Optimization with Piecewise-Linear Decision Rules (RO-PLDR): 
implementing fl24D using piecewise-linear decision rules; 

4. Stochastic Programming (SP): a discretization of fl22p with sampling 
of the uncertainty 5 resulting in a finite number of scenarios. 


In the implementation of the robust optimization models, we considered that 
decision rules only depend on the uncertainty in heat demand and not on 
the realization of the balancing market price. This reduces the modeling of 
the uncertainty to the estimation of price and heat load expectations as well 
as their variance-covariance matrix. Without loss of generality, we made the 
assumption of no autocorrelation of the processes and slightly positive cor¬ 
relation (0.3) between heat load and balancing market prices. The variances 
we used for the heat demand and balancing market prices are consistent with 
the state-of-th e-art on forecasting, which reports RMSE value s of about 7% 
for heat load ( Nielsen and Madsen . 20061: Zelinka et ah . 2013) and 33% for 


balancing market prices fjJonsson et al.l. 120141) . In the case of piecewise linear 


decision rules, we calculated the mean of the positive and negative parts as: 


E{v} = E{r} = !%^ = <Tj 



h 


(29) 


which holds for a zero-mean Normal variable, see Leone et al. f 196ll ). We 
further assume that positive and negative parts of the heat load forecasting 
error contribute equally to the correlation with the balancing market price. 

As far as the stochastic programming model is concerned, we generated 
scenarios by drawing samples from a multivariate Normal distribution fully 
dehned by the characteristics described above. We employed 100 scenarios 
after reduction from an initial set of 2000 perfor med using the fast-forward 
method presented in Heitsch and Romisch f 2003h . We employ scenario fans, 
i.e., a collection of individual paths for the uncertainty branching out at 
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Time period (h) 

Figure 1: Day-ahead electricity price during simulation day in spring 


the root node (representing the hrst-stage) but with no further branches 
at subsequent stages. Note that the non-anticipativity of the solution is 
not fulhlled within this setup. This approximation is needed as the use of 
scenario trees within a multi-period problem with 24 stages is prohibitive. 

The remainder of this section is structured as follows. In Section 14. we 
highlight operational differences between the proposed robust optimization 
model and its deterministic counterpart. Section 14.21 assesses the improve¬ 
ment obtained by switching from linear to piecewise linear decision rules 
within the robust optimization framework. Section |4]3] presents a sensitivity 
analysis showing how the parametrization of the uncertainty set affects the 
performance and the robustness of the solution. Finally in Section 14.41 we 
compare the proposed approach with deterministic optimization and stochas¬ 
tic programming. 

4.1. Daily operation 

Let us consider the operation of the heat and power system during the 
fourth day of the spring week considered in this numerical study. The day- 
ahead electricity price for this day is shown in Figure [H Notably, this price 
peaks during the morning hours. 

The heat production and consumption schedule resulting from the de¬ 
terministic model is illustrated in Figure [2al For each hour, the left bar 
represents supply and the right one demand of heat. Note that the heat 
storage can be on both sides depending on whether it is being discharged 
or charged. Notably, as the heat load is relatively low in this season, the 
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Table 2: Change of dispatch for extraction CHP unit during simulation day 
in spring as a function of T 



Unit 



P 




2 

4 

6 

8 

10 

Total heat dispatch 

MWh 

265.20 

444.27 

622.24 

811.31 

1250.74 

Periods online 

h 

4 

5 

6 

8 

24 


back-pressure CHP unit is sufficient for covering the demand alone (the ex¬ 
traction unit is initially on and quickly turned off). The heat storage is used 
to perform arbitrage hrst by being charged as the back-pressure unit ramps 
up heat and power production in the morning, thus exploiting the electric¬ 
ity price peak, and then by being discharged during the lower power price 
periods. 

The schedule obtained for the same day with the robust optimization 
model with linear decision rules (RO-LDR) is shown in Figure where 
we employed an interval radius equal to 3.2 times the estimated standard 
deviation of the heat demand at each time period and a budget of uncer¬ 
tainty P = 6. Two observations can be drawn from the comparison with 
the deterministic schedule in Figure [2al Firstly, the extra conservativeness 
of the robust optimization approach reduces the extent to which arbitrage 
is performed: notice, for example, the shorter peak of electricity production 
around 11am. The second observation concerns the scheduling of the extrac¬ 
tion CHP unit at the end of the day. Although more expensive than the 
back-pressure unit, this plant is needed to be online for the system to cope 
with possible deviations of the heat demand from the forecast value. 

We conclude this section with Table [21 which shows the total planned 
production from the extraction unit over the 24-hour period as a function 
of the budget of uncertainty P. Notably, as the requirements in terms of 
conservativeness of the solution increase, the higher the schedule for the 
extraction unit. When scheduled, this unit can both ramp up and down in 
response to unexpected deviations of the heat load. 

4 . 2 . Linear vs piecewise linear decision rules 

In the previous section, the robust optimization model with linear decision 
rules is seen to partly replace in the dispatch the back-pressure CHP unit 
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Schedule (MWh) Schedule (MWh) 



(a) Deterministic model 



(b) RO-LDR (interval radius equal to 3.2 times the standard deviation, F = 6) 
Figure 2: Day-ahead dispatch for the system during simulation day in spring 
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Table 3: Change of dispatch for extraction CHP unit during simulation day 
in winter as a function of T for linear and piecewise-linear decision rules 


Decision rule 

Unit 



P 



2 

4 

6 

8 

10 

Linear 

MWh 

7944.00 

7783.77 

7572.89 

7377.16 

7200.06 

Piecewise-linear 

MWh 

7944.00 

7944.00 

7944.00 

7944.00 

7944.00 


with the extraction one, despite the higher cost of the latter, in order to 
guarantee flexibility. 

Table [3] shows that the dispatch decision suggested by the robust op¬ 
timization model with purely linear decision rules can be suboptimal. This 
table illustrates the total heat dispatch for the extraction unit over the fourth 
day of the winter week. Because the peak heat load occurs during this day, 
the system needs both the back-pressure and the extraction CHP plants to 
be online, and must use the peaker to provide the remainder of the heat load. 
Yet, the dispatch of the extraction unit decreases as the uncertainty budget 
P grows. This behavior is caused by the linearity of the decision rule, which 
requires the extraction unit to provide both upward and downward regulation 
in response to heat load deviations. Hence, the setpoint for heat production 
from this unit must be lower than the production capacity if this unit has to 
respond to deviations in heat demand. Piecewise-linear decision rules allow 
for two different re-dispatch strategies depending on whether the heat load 
is lower than expected (in which case the extraction CHP unit reduces pro¬ 
duction) or higher (the peak unit increases production). With purely linear 
decision rules, these two re-dispatch strategies must be one and the same (the 
extraction CHP unit is used to cover both positive and negative deviations 
of the heat load). 

The suboptimality of the dispatch described above has a signihcant im¬ 
pact in terms of proht for the system. Figure [3a] shows the improvement in the 
objective function value obtained when switching from linear to piecewise- 
linear decision rules. Notably, this theoretical improvement is rather large 
and increases with the required conservativeness of the solution, both in 
terms of size of the interval and budget P for the uncertainty set dehned in 
Section 1221 The improvement just mentioned, though, is only theoretical as 
it is calculated on the assumption that the decision-maker will follow linear or 
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Figure 3: Improvement in expected profit when switching from linear to 
piecewise-linear decision rules on winter day 


piecewise-linear decision rules in practice, which is not necessarily true. An 
approximation of the practical improvement can be determined by hxing the 
unit commitment and dispatch solution of the robust optimization model in 
a stochastic programming model implementing a discretized version of fl22p . 
The solution is then evaluated out-of-sample by considering the redispatch 
(not necessarily affine on the heat load deviation) over a set of 100 scenarios. 
Average results for the same day are shown in Figure [3bl Notably, the im¬ 
provement is lower but still signihcant. However, it should be remarked that 
the stochastic programming model does not preserve the non-anticipativity 
of the solution, as discussed at the beginning of this section. Hence, this also 
represents an approximation of the actual improvement of piecewise-linear 
decision rules. 

4-3. Assessing the priee of robustness 

One aspect of critical importance in robust optimization is hne-tuning 
the parameters governing the size of the uncertainty set. A certain level 
of robustness is required by the very nature of the problem at hand. In¬ 
deed, failure to meet heat demand may have some important consequences, 
either hnancial or in terms of corporate image. On the other hand, increas- 
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(a) Average profit (b) Heat-load not served (worst case) 

Figure 4: Out-of-sample average profit and heat-load not served obtained 
with the robust optimization model under different parameters on winter 
day 


ing robustness requirements are met at the expense of the overall financial 
performance (in expectation). 

Figure |4] illustrates the trade-off between average proht and the worst- 
case realization of heat-load not served during the winter day characterized 
by peak heat load. The results in the hgure are obtained from out-of-sample 
simulation where the unit commitment and dispatch solution is fixed in the 
stochastic programming model. Notably, the average profit decreases as the 
conservativeness requirement (i.e., interval size and uncertainty budget) is 
strengthened in Figure HU However, increasing the robustness requirement is 
beneficial only up to a certain point in terms of heat-load not served. Indeed 
as Figure Hb] confirms, there appears to be no need to push the interval radius 
of the uncertainty set much beyond three times the standard deviation. 

For the sake of brevity, we omit similar plots obtained for spring, sum¬ 
mer and autumn, where we observe similar trends with in general lower 
requirements in terms of conservativeness. Note that those results are partly 
included in the following section. 

Robust optimization vs alternative methods 

This section is dedicated to a comparison of the robust optimization ap¬ 
proach with two other state-of-the-art techniques described at the beginning 
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of this section: deterministic optimization and stochastic programming. We 
evaluate the unit commitment and dispatch solutions from the different mod¬ 
els by hxing them in a stochastic programming approach and evaluating their 
results after re-sampling the uncertainty (i.e., with different scenarios than 
the ones used in the optimization of the stochastic programming model). 

Results in terms of average proht and heat-load not served (abbreviated 
as “heat LNS”) from the winter, spring, summer and autumn week are shown 
in Tables m |5l |6]and[71 respectively. For the robust optimization approach we 
only included parameter choices that result in Pareto-optimal points in the 
space of expected proht and worst-case heat-load not served. The average 
proht does not include penalties for shedding heat load, nor it discounts the 
fuel saved by the operator of the system as a result of curtailments. 

Both the robust optimization and the stochastic programming approach, 
as expected, signihcantly reduce the heat-load not served obtained with the 
deterministic optimization model. Furthermore, if the average heat-load not 
served is discounted from the average proht with the sum of avoided fuel cost 
and a penalty (or the fuel cost of an expensive backup unit), both models 
outperform the deterministic approach. This is particularly true in winter 
and spring, where the expected heat load not served is rather high for the 
deterministic solution, see Tables |4] and [5l 

Moreover, the robust optimization approach with piecewise-linear deci¬ 
sion rules provides a hexible framework where conservativeness of unit com¬ 
mitment and dispatch can easily be traded-oh with average hnancial perfor¬ 
mance simply by tuning the interval size and the budget of the uncertainty 
set. For sufficiently large values of these parameters, heat-load shedding can 
be completely removed in all the cases considered. This is particularly im¬ 
portant in comparison to the stochastic programming approach, which incurs 
some load shedding in all weeks excluding the summer one. We remark that 
this heat-load shedding is caused by the necessary atomization of the sample 
space of the uncertainty, despite the use of a large sample set with a scenario 
reduction technique. We deem this drawback rather critical, as it would not 
be possible to get rid of the residual load shedding by considering a risk 
measure in the stochastic programming formulation. 

The results in Tables HHT] show that the stochastic programming solution 
achieves higher average prohts for comparable levels of heat-load not served. 
This behavior can be explained at least partly by the fact that stochastic pro¬ 
gramming is used to evaluate all the models. Indeed, this evaluation scheme 
does not reward the non-anticipativity of the solution, which robust opti- 
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Table 4: Out-of-sample simulation results on winter day. Only Pareto- 
optimal versions of the robust optimization approach are shown 


Model specihcation 


Avg. proht (€) 

Heat LNS (MWh) 

method 

radius {4^ sd) 

r(~) 

largest 

expected 


2.00 

2 

449 301.37 

329.38 

7.88 


2.00 

4 

446 337.45 

198.93 

3.38 


2.00 

6 

443 598.68 

144.68 

2.51 

RO-PLDR 

2.40 

2 

443 474.44 

127.20 

0.51 

2.40 

4 

440 776.45 

39.64 

0.29 


2.40 

8, 10 

438 086.22 

31.58 

0.20 


2.80 

2 

432 707.89 

16.59 

0.05 


3.20 

2-10 

419 281.79 

0.00 

0.00 

DET 

- 

- 

455 815.78 

745.78 

69.52 

SP 

- 

- 

448 945.73 

159.90 

2.13 


mization models with decision rules guarantee (on the contrary of stochastic 
programming), and thus it penalizes the former. 

5. Conclusions 

In this paper, we address the determination of the optimal day-ahead 
unit commitment and dispatch of systems including heat and power produc¬ 
tion units as well as heat storages. Owners of such units face this problem 
as they have to decide their trading strategy on the day-ahead power mar¬ 
ket and communicate their schedules for heat production in advance. As a 
result of the look-ahead time, this optimization problem is subject to the 
uncertainty in heat demand and in electricity prices, as these are unknown 
at the time of decision-making. Furthermore, the problem is inherently dy¬ 
namic (and multi-stage) because of the intertemporal constraints governing 
the production and storage units. 

We propose a formulation of the problem based on robust optimization 
where optimal recourse decisions are approximated by linear or piecewise- 
linear functions of the uncertain parameters. Uncertainty is accounted for 
in this formulation through the dehnition of an uncertainty set as well as 
of (conditional) means and variance-covariance matrices for the uncertain 
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Table 5: Out-of-sample simulation results on spring day. Only Pareto- 
optimal versions of the robust optimization approach are shown 


Model specihcation 


Avg. proht (€) 

Heat LNS (MWh) 

method 

radius (7^ sd) 

TB 

largest 

expected 


2.00 

2 

251514.22 

217.26 

5.70 

RO-PLDR 

2.40 

2 

250 505.55 

6.78 

0.01 


4.00 

2 

248484.94 

0.00 

0.00 

DET 

- 

- 

260 927.88 

340.15 

49.74 

SP 

— 

— 

257295.35 

76.17 

1.54 


Table 6: Out-of-sample simulation results on summer day. Only Pareto- 
optimal versions of the robust optimization approach are shown 

Model specihcation , n, 

Avp-. profit 

Heat LNS (MWh) 

method radius (# sd) P (-) 

largest expected 

RO-PLDR 2.00 2 88423.27 

DET - - 94413.41 

SP - - 92 356.87 

0.00 0.00 

117.00 16.10 

0.00 0.00 


Table 7: Out-of-sample simulation results on autumn day. Only Pareto- 
optimal versions of the robust optimization approach are shown 

Model specihcation 

Avp- proht (^) 

Heat LNS (MWh) 

method radius (7)^ sd) 

P (-) 

largest 

expected 

RO-PLDR 

2.40 

DET 

SP 

2.00 249 631.81 

2.00 248 673.51 

259135.23 

255 836.17 

90.35 

0.00 

238.72 

99.71 

0.41 

0.00 

12.01 

0.54 
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parameters, i.e., heat demand and power prices. Furthermore, affine decision 
rules allow the modeler to enforce dependence of each recourse decision solely 
on the subset of stochastic parameters whose actual value unfolds before the 
implementation of the recourse decision itself. In other words, the proposed 
model complies with the non-anticipativity condition in multi-stage decision¬ 
making processes. 

An extensive numerical study based on data from the Copenhagen area 
in Denmark highlights several properties of the proposed model. Through an 
example of daily operation, we show that the higher requirements in terms 
of conservativeness of the solution (i.e., its capability to cope with unfore¬ 
seen heat-load deviations) in the proposed model slightly reduce the pos¬ 
sibility of exploiting the flexibility of the heat system to perform arbitrage 
in the electricity market. Further financial results from the study suggest 
that a significant improvement can be obtained by switching from linear to 
piecewise-linear decision rules. Moreover, we perform a sensitivity study as¬ 
sessing the trade-off between average profit and conservativeness in terms of 
heat-load not served resulting from tuning the parameters (size and budget) 
defining the uncertainty set. Finally, we compare the proposed model with 
the alternative approaches of deterministic optimization and stochastic pro¬ 
gramming. We show that, differently from the competing approaches, our 
model (if appropriately tuned) can result in complete immunization from 
the heat-load uncertainty. Under a mild penalization of heat-load not served 
(e.g., including the activation cost of backup units), the average profit can 
improve that obtained with deterministic optimization. If solely consider¬ 
ing average profit, stochastic programming seems to outperform our model. 
However, this evaluation is performed via stochastic programming (after re¬ 
sampling the uncertain parameters from the assumed distribution), which 
does not enforce non-anticipativity of the recourse decisions. Devising (and 
testing the models on) a non-anticipative evaluation scheme is an interesting 
direction for future research. 

Other interesting topics for future research include considering piecewise- 
linear decision rules with more than two pieces, the inclusion of binary vari¬ 
ables in the recourse decisions, the use of more sophisticated models of the 
uncertainty (e.g., based on time-series models) and the inclusion of more 
market stages in the model (e.g., reserve or intra-day electricity markets). 
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